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ABSTRACT 

The formation and collapse of a protostar involves the simultaneous infall and outflow 
of material in the presence of magnetic fields, self-gravity, and rotation. We use self- 
similar techniques to self-consistently model the anisotropic collapse and outflow by a 
set of angle-separated self-similar equations. The outflow is quite strong in our model, 
with the velocity increasing in proportion to radius, and material formally escaping 
to infinity in the finite time required for the central singularity to develop. 

Analytically tractable collapse models have been limited mainly to spherically 
symmetric collapse, with neither magnetic field nor rotation. Other analyses usually 
employ extensive numerical simulations, or either perturbative or quasistatic tech- 
niques. Our model is unique as an exact solution to the non-stationary equations of 
self-gravitating MHD, which features co-existing regions of infall and outflow. 

The velocity and magnetic topology of our model is quadrupolar, although dipolar 
solutions may also exist. We provide a qualitative model for the origin and subsequent 
evolution of such a state. However, a central singularity forms at late times, and 
we expect the late time behaviour to be dominated by the singularity rather than 
to depend on the details of its initial state. Our solution may, therefore, have the 
character of an attractor among a much more general class of self-similarity. 
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1 INTRODUCTION 

It is clear that outflows often co-exist with infall as proto- 
stars form within the collapsing cores of molecular clouds 
(Bertout 1989; Andre et al. 1993). Infall and outflow both 
appear to be present for much of the protostellar "main se- 
quence," from rapidly accreting embedded Class objects 
to fully formed T Tauri stars. This suggests that the dy- 
namics leading to the formation of a protostar are more 
complex than simple radial infall, and are dominated by 
strongly anisotropic motions. We present a new model for 
the anisotropic collapse of a molecular cloud core, which self- 
consistently treats the effects due to self-gravity, magnetic 
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fields, and rotation, as the central protostellar core grows 
and a bipolar outflow develops. 

Self-gravitating models of protostellar collapse have 
usually been limited to the classical solutions with spher- 
ical symmetry (Larson, 1969; Penston, 1969; Shu, 1977), in- 
cluding the elaborations and clarifications in related works 
(Hunter, 1977; Whitworth and Summers, 1985; Henriksen, 
Andre and Bontemps, 1997, hereafter HAB1997). GaUi and 
Shu (1993, hereafter GS) presented a very interesting calcu- 
lation, which included the effects of anisotropy as a per- 
turbation about the classical Shu (1977) inside-out col- 
lapse solution. Subsequently, Li and Shu (1996) presented 
a quasi-static calculation. An approximate analytic self- 
similar treatment based on a dynamic termination of the 
ambipolar diffusion models has also been given recently and 
developed to the point of comparison with observations (e.g. 
Basu, 1997). However the bipolar outflow was not integral 
to any of these papers, as it is in the case of the present 
work. 

We have previously studied steady-state solutions for 
simultaneous infall and outflow late in the evolutionary se- 
quence, after the dominant central mass had already formed 
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(Henriksen and Valls-Gabaud, 1994; Fiege & Henriksen 
1996, hereafter FHl; Lery, Henriksen, & Fiege 1999, here- 
after LHF). Thus, our previous models apply only to very 
late times in the formation of a protostar. The philosophy 
of these articles was that bipolarity could be studied using 
scale free solutions near the central singularity without wor- 
rying too much about the initial state of the flow. This is 
similar in spirit to the development of the Larson-Penston 
self-similar solution from non-self-similar initial conditions. 

The model presented in this paper takes quite a dif- 
ferent approach by treating the time-dependent problem of 
accretion and simultaneous outflow in a dynamically collaps- 
ing and self-gravitating core. Thus, we study an early stage 
of stellar formation when the star has not yet formed, and 
most of the gas still resides in the surroundings. A limit to 
the self-gravitating regime is certainly set when the mass of 
the central object dominates that of the surroundings. Thus 
our present solution is a natural complement to our earlier 
studies. 

The present model is best described as an inner "set- 
tling" solution, which follows the assembly of the protostar 
in detail. It is limited to smaller spatial scales than our pre- 
vious work, but the flow structure that we predict would 
presumably be embedded within a larger collapsing region. 
This larger region might include nearly steady-state outflows 
and an accretion region of the self-similar type that we have 
previously discussed in FHl and LHF (See Section ^). 

Numerical simulations have also been used to 
study non-isotropic self-gravitating collapse. For example, 
Tomisaka (1998) used a multi-grid MHD code to study the 
gravitational collapse of a segment of a magnetized and 
slowly rotating filament. His simulation resulted in a rotat- 
ing pseudo-disc, which produced an outflow as the central 
object grew. Most interesting, from our point of view, is 
that the late stages of the calculation were dominated by 
a quadrupolar velocity field. Such quadrupolarity can de- 
velop in super- Alfvenic flows when infalling material near 
the midplane is deflected up the axis due to a combination 
of pressure gradients, magnetic fields, and the centrifugal 
barrier. Our present calculation shows that this mechanism 
can operate on smaller scales between the growing boundary 
of the hydrostatic core and the pressure dominated region 
external to the region of self-similarity. The development of 
quadrupolar structure and the connection of our model to 
the exterior region is discussed in Section |5[ 

Our central assumption is that this early stage of star 
formation is dynamic rather than quasi-static. This is sug- 
gested by various lines of evidence (Basu, 1997; Foster and 
Chevalier, 1993; HAB, 1997), although this assumption can- 
not be regarded as certain (Basu and Mouschovias 1995a; 
Basu and Mouschovias 1995b; Galli and Shu, 1993; Li and 
Shu, 1997). In any case, the central regions must ultimately 
become hydrostatic to allow for the growth of the naissant 
protostar. Thus, we expect and do indeed find "settling" 
solutions, in which the radial velocity goes to zero near the 
centre, where the infalling material actually forms the stellar 
core. 

In the following section, we derive our basic self-similar 
equations from self-gravitating MHD, under the assumption 
of self-similar flow. We show, in Section ^ that our equations 
admit an exact and completely analytic class of solutions. 
We use these solutions to derive several interesting analytic 



results, which illustrate and constrain the properties of the 
model. Our most important result is that substantial out- 
flow velocities can be obtained without resorting to heating 
the material, as in FHl and LHF. The axial outflow veloc- 
ity is never more than twice the equatorial inflow velocity 
on a sphere of some given radius, at some instant of time. 
However, we find that the outflow velocity increases linearly 
with distance, so that substantial velocities are obtained far 
from the origin. By following the motion of individual fluid 
elements, we find that all such fluid elements escape to ra- 
dial inflnity in the finite time required for the central sin- 
gularity to develop. This escaping gas would presumably in- 
teract with the external medium in a complicated manner, 
which we discuss in Section Our Discussion section also 
presents a simple, qualitative model to provide one possi- 
ble way in which the quadrupolar structure might arise. We 
also note that models with dipolar geometry are also possi- 
ble, although we do not find them explicitly in the present 
calculation. 

Finally, we note that a preliminary account of this work 
was presented at the Cracow meeting on "Plasma Turbu- 
lence and Energetic Particles in Astrophysics" (1999, M. Os- 
trowski and R. Schlickeiser, eds.). 



2 SELF-SIMILAR ACCRETION AND 
OUTFLOW 

We assume that the protostellar development begins with 
the primarily radial collapse of a cold molecular core. In re- 
ality, cores are not really expected to be spherical (Tomisaka 
et al. 1988, Myers et al. 1991, Ryden 1996, Fiege & Pudritz 
2000a), but all that is really required is that there be a 
substantial collapse in a plane perpendicular to the initial 
magnetic and rotation axes (assumed to be parallel for sim- 
plicity). The core may have been in equilibrium initially as 
a singular isothermal sphere (SIS) as in GS, or it may have 
had a more complicated internal structure as suggested by 
HAB. The initial conditions provide a set of characteristic 
scales, which suggest the set of conserved quantities used to 
define the class of self-similar symmetry. For example, the 
obvious constants for an unbounded SIS are the sound speed 
Cs and Newton's constant G. GS used these constants to de- 
fine a "class" of self-similarity obeyed during the collapse. 
On the other hand, a constant external pressure bounding 
a truncated self-gravitating sphere, together with G, would 
define a different class of self- similarity. Yet another similar- 
ity class would be appropriate if there were a characteristic 
time due to rotation. 

Carter and Henriksen (1991) developed a mathemat- 
ical formalism for determining the most general class of 
self-similarity possible, which naturally includes all possi- 
ble initial states. This method has been used successfully in 
studying the evolution of coUisionless n-body systems (see 
e.g. Henriksen, 1997), but is equally well-suited for study- 
ing the dynamical collapse of a magnetized, self-gravitating 
core. This technique can be used to demonstrate the exis- 
tence of the separable (in radius, poloidal angle, and time) 
"settling" solution presented here, as a special case within a 
more general class of self-similarity. It is possible that these 
separable solution might represent an "attractor" within this 
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larger class of self-similar models (Aburihan, 1999), but we 
do not prove this here. 

It is important to address the question of how the 
quadrupolar magnetic and flow structure assumed by our 
model might originate. Such questions are really beyond 
the scope of our self-similar treatment, since self-similar so- 
lutions are often intermediate, in the sense that they are 
disconnected from their boundary conditions in space or 
time (see discussion in FHl). Thus, our ideas regarding the 
origin of our self-similar model are somewhat speculative. 
Nevertheless, we suggest a reasonable scenario in Section |[ 
by which a rotating and magnetized cloud could evolve a 
quadrupolar flow and magnetic field structure during col- 
lapse. This would arise as a consequence of poloidal pres- 
sure gradients, the centrifugal barrier encountered by the 
collapsing cloud, and magnetic reconnection effects which 
change the topology of the field in regions where localized 
field reversals occur. We discuss the details of this scenario 
in Section |5|, and turn now to a derivation of the basic equa- 
tions. 



2.1 Equations 

The basic non-dimensional quantities from which we con- 
struct our self-similar model are given by the poloidal angle 
6 and the variable 



X 



r 

Cs t' 



(1) 



where r is spherical radius, t is time, and Cs is a fiducial 
sound speed. The local sound speed need not be constant; 
the constant Ca in equation ^ is only meant to be typical of 
the initial conditions. The minus sign is included to make 
X positive definite, since our model starts with t large and 
negative, evolving toward a singularity at t = 0. Note that 
our X is essentially the same as the self-similar variable 
used by Shu (1977) and more recently by GS in the con- 
text of collapsing isothermal spheres, although their model 
is strictly isothermal. We provide a self-similar framework 
in this section which would, in principle, allow one to follow 
the collapse of a SIS. 

Our treatment admits considerable freedom in choos- 
ing the equation of state (EOS). One particularly simple 
non-isothermal choice, which we discuss later in this section 
and use extensively throughout this paper, makes the par- 
tial differential equations separable, which allows us to find a 
singular and completely analytic solution to our self-similar 
equations. Unfortunately, this particular solution does not 
match on to the SIS at early times. A more complete treat- 
ment of the collapse problem would solve (with greater ef- 
fort) the general self-similar PDEs directly, beginning with 
realistic initial conditions, such as a SIS threaded by a mag- 
netic field, rather than seeking out special separable forms. 
However, our singular solution might represent an "attrac- 
tor" among the more general class of self-similarity, which 
would make it a valid end state for a wide variety of initial 
conditions. This possibility is further explored in Section ^ 

The actual temperature and sound speed in our model 
vaxy as functions of r, 8, and i in a way that is uniquely de- 
termined by the self-similarity. Generally speaking, the gas 
heats up during collapse in our model, although with sig- 
nificant temperature structure in 9. Isothermality in cores 



is maintained primarily by the efficient cooling provided 
by molecular lines. An anisotropically collapsing core would 
likely evolve toward a more complex temperature structure 
once the dynamical timescale r/ca becomes shorter than the 
cooling time. Thus, we expect the self-similarity to develop 
from the inside out, with the central regions of the collapsing 
core evolving more rapidly toward a self-similar flow pattern, 
while the outer layers remain nearly isothermal until later 
in the collapse. 

The most general class of self-similarity based on X and 
6 is expressed by writing all physical quantities as functions 
of these variables: 



B 

$ 
P 



CsW{X,e) 
^b(X,6l) 



G r 



c; I po{t) pix, 

47rG ' " ^ 



(2) 
(3) 

(4) 
(5) 
(6) 



The notation for the dimensional magnetohydrodynamic 
(MHD) quantities on the left land side of these equations 
is standard, and the quantities p,, W, b, if>, and p on the 
right hand side are dimensionless forms of the density, veloc- 
ity, magnetic field, gravitational potential, and "dynamical" 
pressure respectively. Note that we separate the pressure 
into a "dynamical" pressure term p, and a background pres- 
sure term po(t). The term involving po(t) has no dynamical 
consequences whatsoever, since it does not contribute to the 
pressure gradient. It is necessary because realistic solutions 
are obtained only when the "dynamical" part of t he pres- 
sure p < 0, as we shall further discuss in Section 4.1. Finally, 



note that the overbars distinguish these quantities from the 
separated forms of these variables, presented later in this 
section. 

We now use the preceding ansatz of self-similarity in 
each of the equations of self-gravitating MHD to obtain our 
general set of self-similar equations. Poisson's equation, the 
continuity equation, and the condition that there are no 
magnetic monopoles are written respectively as follows: 
Poisson's Equation 



Continuity 



(7) 



(Wr + X)Xdxfi + fiXdxWr + ^doifiWo sin 6) = 0; (8) 
No Magnetic Poles 



dxiXbr) + -—deismObe) = 0. 



(9) 



The 3 components of the MHD induction equation are given 
by the following: 
Induction Equation - f 



X^dxbr 



-dg{sm OE^) 



Induction Equation - 6 

- Xdxbe + dxEs = ; 



(10) 



(11) 



nnnn o A O TV/fMO A C r\r\r\ 



4 Mahmoud Ahurihan, J.D. Fiege, R.N. Henriksen & T. Levy 



Induction Equation - (j) 

- Xdxb^ = dxEg - -^deEr 



(12) 



where Er, Eg, and E^ are related to the electric field and 
given by 



E^ = ~{Wrbg - Wgbr), 

Eg = -{W4,br-Wrh), 

Er = ^{Wgb^-W^bg). 



(13) 
(14) 
(15) 



Finally, the 3 components of the momentum equation are 
given as follows: 
Momentum - f 

X{Wr + X)dxWr + WedeWr ~ (Wi + Wl) 

^ ^ _ }-XdxP - Xdx^ + -bedebr 
MM M 

-^Xdx{bl+bl) ■ 

Momentum - 6 

X{Wr + X)dxWg + WgdgWe + WrWg - cot ( 
1_ -1- - 1-2 

= ——dep — deip + —brXdxbe — —bs cot 9 
M M M 



(16) 



1 

'2fi 



(17) 



Momentum - 



1 



X(Wr + X)dxW^ + —Wgdg{smeW^) + W^Wr 

= -WXdxb^ + ^-6eae(sine6^). (18) 
fi sm o 

Note that the self-similar forms presented in equations 
^ to 1^ would represent a SIS 'd fi — 2 and W = 0. A complete 
treatment of the self-similar problem would involve solv- 
ing the {X,9) partial differential equations presented above, 
starting from a slightly perturbed SIS at some initial time to 
{to = — oo for an infinitestimal perturbation). In principle, 
this approach would allow us to follow the entire anisotropic 
development of the SIS all the way to the development of 
the singularity at t = 0. This could be done under various 
assumptions regarding the initial conditions for the rotation 
and magnetic field, resulting in a very complete model of 
protostellar collapse. This ambitious project will be left for 
future work. Instead, we take a more modest approach in 
this paper, by solving the restricted problem where all of 
the variables fi written with overbars in equations ^ to ^ 
can be written in the separable form 



Mx,e) = hx'''Me), 



(19) 



where the ki are appropriate normalizing constants. In the 
case of p, and p these constants are each Att,. They are equal 
to -|-1 elsewhere. The factor that depends on X is eliminated 
by carefully balancing powers of a^, resulting in the following 
system of ordinary differential equations in 6, written here 
in the same order as in equations |^ to |l^ and |l6| to |l^. 



6* + ^—rdg(sm9dg^) = 4nu,, 
sm 

SWr + 2+ J—dg (sin eWg) + Wgdg \n^I^O, 
smd 



(20) 
(21) 



36, + -—dg{sm 6be) = 0, (22) 
sm 

- 2br + -^de [sme{Wrbe - Wgbr)] = 0, (23) 

2bg - 3{Wgbr - Wrbg) = 0, (24) 

- 2b^ + 3{W^br - Wrb^) + dg{W4,bg - Wgb^) = 0, (25) 



Wr + Wr + WgdgWr " (W! + W^) 

-2P-2^+-^\{bgdgbr~2{bj+bl)], (26) 

H 'inn 



Wg + WgdgWg + 2WrWg -WlcotO = ~ de* 



47r/i 



2brbg - -dg{br + 6^) - b^ cot f 



Wg 



W4, + 2WrW^ + ^dg{smeW4,) 



1 



2brbj, + -^—-dg{sm9bs) 



(27) 



(28) 



Equations ^ to ^ are the final equations, whose solutions 
we study for the remainder of this paper. 

The physical variables are given explicity by the follow- 
ing forms, in which we have replaced X by r and t, using 
the definition of X given in equation |l[ 



V 

B 

P 
P 
$ 



^/Gt^ 
1 

1 



w(e), 
b(e). 



rhi + r^pie)] 



t2 



*(e). 



(29) 
(30) 

(31) 
(32) 

(33) 



Note that we have replaced the term involving po(i) in equa- 
tion ^ with one involving only the constant po, which can 
be chosen a posteriori to keep the total pressure positive 
throughout the region of interest. Formally we have set 
Po{t) = ^rtp*o. 

One should realize that our final set of equations can 
also be obtained by using these physical forms directly in 
the equations of self-gravitating MHD, without ever writ- 
ing down the {X, 9) forms of the equations (equations |^ to 
p^). However, our treatment has the advantage of provid- 
ing the equations within the context of a broader class of 
self-similarity, which may be useful for future analysis. It 
is important to realize that Cs has disappeared from these 
equations, as would any dimensional constant used in equa- 
tion ^. Thus there is a "route" to these equations from quite 
general initial conditions, which is the main reason for our 
suspicion that it may be an attractor. 

The nature of our model as an internal settling solu- 
tion is clear from the self-similar forms in equations ^ to 
before we even solve the equations. We observe that all 
velocity components go to zero proportionally to r, which 
also implies that there is rigid rotation on each cone defined 
hj 6 = const. Of course, the angular velocity is allowed to 
vary between cones, and material flows between them. 
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A difficulty arises when one tries to solve the self-similar 
form of Poisson's equation given in equation Any solution 
to equation po| must satisfy the boundary conditions ^''(0) = 
and »I''(7r/2) = 0, so that there is no 9 component of the 
gravitational acceleration at either boundary. The solution 
to the homogeneous form of the equation (jj, — 0) which 
satisfies these boundary conditions, is the following: 



*/iom = cPa (cos e) = c- (3 cos^ 6-1) 



(34) 



where P2{cos8) is the second order Legendre polynomial 
and c is an arbitrary constant. Note that the homogeneous 
solution, with arbitrary c may be added to any particular 
solution of equation |20| which produces an additional term 
in the gravitational acceleration. The solution to Poisson's 
equation is non-unique for our system, since any such so- 
lution satisfies the boundary conditions at the polar and 
equatorial boundaries. This non- uniqueness originates from 
the self-similarity. Poisson's equation has a unique solution 
provided that appropriate boundary conditions can be spec- 
ified everywhere on a surface enclosing the region of interest. 
However, there is no unique way to specify the potential at 
radial infinity for our distribution of matter, which in prin- 
ciple extends to infinity. In practice our solution must be 
halted at some surface R{9, t) whose form depends on the 
matching to an external medium. The best that we can do 
is to limit ourselves to the special case where ^(6) is con- 
stant and the density distribution is spherically symmetric 
at all times. In that case, ^' should be contant on physical 
grounds, so that the isopotential surfaces are spherical as 
well. Poisson's equation (^0|) is then trivially solved: 

* = I^M- (35) 

We use this equation and the underlying assumptions of 
constant /i and ^ for the remainder of the analysis in this 
paper. 

Note that the restriction that ^ = const implies a sort 
of "poloidal incompressibility" on our system, which takes 
the place of a thermodynamic equation of state (EOS) in our 
system, in the same way that true incompressibility replaces 
a thermodynamic EOS in incompressible fluid dynamics. It 
is generally not possible to impose any additional EOS di- 
rectly relating p(9) to ^(9) or any other variable, without 
imposing complete boundary conditions on Poisson's equa- 
tion. 

We commented in Section |^ that our solution must only 
be valid near the centre of a collapsing core, where the dy- 
namical timescale is much shorter than the cooling time of 
the molecular gas. Realistically, our settling solution must be 
embedded within a nearly isothermal exterior region, which 
would remain nearly isothermal, as a result of the efficient 
cooling provided by molecular lines and the longer dynami- 
cal timescale outside of the collapsing region. The boundary 
joining our model to such a region would undoubtedly be 
complex in both shape and internal structure, with shocks 
arising near the outfiow axis. A more complete model, which 
explicitly includes this external region would be very difficult 
to treat analyticaUy. The equatorial region is expected to 
contain both acretion discs and magnetically neutral points 
as suggested above, while there may be more violent activ- 
ity (including shocks in super Alfvenic fiow) near the axes. 
None of this can appear in the simple asymptotic forms that 



we study here. The only effect of the external region in the 
present calculation is through the background pressure po (t) 
in equation which has no effect on the dynamics what- 
soever. 

The self-similar forms given by equations ^ to ^ be- 
come singular a.t t = Q over the entire domain of validity, 
and they cannot be continued beyond this singularity. This 
is unlike the usual point singularity (i.e. a "core" of finite 
mass) that forms at t = r = in the spherically symmetric 
collapse models, across which the solution may be contin- 
ued into the core accretion phase. Such a solution (but with 
appropriate asymmetry) would be external to the solution 
presented here. 

The inner boundary of our solution is the growing hy- 
drostatic core and realistically, we expect the self-similar so- 
lution to vanish altogether before t = Q. This can happen if 
the outer boundary is shrinking in time as the inner bound- 
ary grows. We can make this plausible near the equatorial 
plane by observing that outside the transition region we can 
expect a mass flux cx c^/G, as given by Shu (1977). At the 
"boundary" R of our inner region equations (^ and (^ show 
that the mass flux scales as /t^ and thus by equating the 
two expressions we obtain R oc Cst. Consequently the outer 
boundary of our solution may be expected to shrink onto the 
hydrostatic core before t ^ 0, thus removing the mathemat- 
ical singularity from the domain over which our solution is 
valid. 

We turn in the next section to an analysis of these equa- 
tions. 



3 INTEGRALS 

There is considerable redundancy in our equations. We note 
that equation 02 can be derived from equations §1 and 1^ 
by simple algebra. There are also several integrals that can 
be derived from our equations. Equation ^is already in the 
form of an integral. Combined with equations ^ and we 
easily find the following two integrals, in simplified form: 



be 

br 



= qWe, 
= q(Wr 



2 

+ 3 



(36) 
(37) 



where g is a constant. A third integral can be obtained 
by inserting the first two into equations]^ to and seeking 
a relation between 6^ and Wcj,: 



(38) 



where Q is another constant. Note that these integrals allow 
us to remove the self-similar magnetic field b entirely from 
the equations to be solved. 

Although it is not essential to the arguments of the 
present paper we might reflect a little on the general signif- 
icance of these integrals. They are likely to be more general 
in fact than our particular solution. In their present physical 
form they read: 

/2 r 



B 



VGt \St 



Be = — 



Bif, — 



VGt 
Qr sin 9 



Gt^ 



VGt 



Vtf,. 
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If we rearrange these equations into vector form as 



2 r 
V 

3 t 



B + - xr, 

3 qt 



where f2 is along the axis of symmetry, then we can recognize 
a kind of Ferraro's theorem (in the form v = const, x B+lj x r 
where uj is the angular velocity of the field line) with time 
dependence. The constants are time-dependent here and the 
material velocity is relative to a freely-falling, zero-energy, 
Keplerian observer. One might expect the time dependence 
to be different for different self-similar symmetries. 



4 ANALYTIC SOLUTION AND ANALYSIS 

4.1 Boundary Conditions and Method of Solution 

We now turn to the task of actually solving the self-similar 
equations presented in the previous section, subject to all 
appropriate boundary conditions. With the help of the in- 
tegrals given in equations |36| to ^ equations ^ to ^ can 
be reduced to a system of equations involving only the self- 
similar velocity components W, pressure p, and their deriva- 
tives with respect to 6. Solving for the derivatives, the equa- 
tions can be written as a dynamical system in standard form, 
which specifies a boundary value problem in four variables 
on the interval ^e[0, vr/2]. 

The boundary conditions at the symmetry axis and 
equatorial plane are as follows. The self-similar velocity com- 
ponents We and must vanish at ^ = 0, where we also 
require that Wr > for an outfiow solution. The velocity 
component We must vanish at the equatorial plane 6* = 7r/2 
for quadrupolar symmetry. The "dynamical" component of 
the pressure p must be negative everywhere, since the pres- 
sure would increase radially outward otherwise, according 
to the self-similar form for the pressure given by equation 
The total pressure can be made positive throughout the 
region of interest, which is bounded on the inside by the 
forming hydrostatic core and on the outside by the external 
molecular gas that provides the source of the background 
pressure poit) in our model. Note that we have demanded 
strict analyticity everywhere on the angular domain for the 
solutions presented in this paper. This restriction could, in 
principle, be relaxed by allowing a singular radial velocity 
along the axis, in the spirit of FHl and LHF. 

The following symmetry conditions allow us to extend 
our solutions from [0,7r/2] to the full sphere [0,7r]. Under 
reflection about the midplane (9 n — 6), we assume that 
We — > —We, hr — > — 6r, and — > ~h^. All other quantities 
remain unchanged upon reflection. Note that the reversal of 
the magnetic field at the midplane requires the existence of a 
current sheet in the equatorial plane, which is in accordance 
with the scenario discussed in Section |5| for the generation 
of the quadrupolar field. Note that the symmetry assumed 
here differs from our previous work, where we assumed that 
We — » —We and be — > —be, with all other quantities un- 
changed upon reflection. In this case, continuity requires 
that both We and be must vanish on the midplane, thus 
requiring quadrupolar symmetry for all solutions. This is 
not required in general for the present model, which admits 
solutions with dipolar symmetry in principle. Different equa- 
torial boundary conditions are required, notably be and 



continuous, but this is permitted in principle by the time de- 
pendence. We will, however, concentrate on the quadrupolar 
class of solutions in this paper. 

We proceeded flrst with a numerical survey of the solu- 
tions in an extended region of parameter space. However, it 
gradually became clear from the extremely simple appear- 
ance of the solutions and one extremely robust result (see 
equation |5^ below), that an underlying analytical solution 
was at play. Once we realized this, it was relatively easy to 
discover the form of the analytical solution. Given the pa- 
rameters ^, q, and il, we assume a solution of the following 
form and solve for the constants ci through cg in the result- 
ing algebraic set of equations, which are derived from the 
reduced system of equations involving only the self-similar 
velocity components and pressure: 



Wrie) 
Weie) 

P{0) 



Cl -|- C2 cos 

C3 sin(26') 

C4 sin 6 

Co + ca sin^ 



(39) 
(40) 
(41) 
(42) 



After straightforward but tedious algebraic manipulations, 
it turns out that the constants C2 through ce can all be 
expressed in terms of ci and the other three parameters 
mentioned above: 



C2 
C3 

C4 

C5 
C6 



-2 - 3ci 

£2 

' 2 



C2qQ 



6Tvfi{l + 2ci) +q^C2 
-7rAt^-/i(l + 2ci)(l+ci) 

■!- {-{C4q + + 2nfi [cl - C2(l + Ci)] } 



47r 



(43) 
(44) 

(45) 

(46) 

(47) 



The parameter space is four-dimensional, with solutions to 
our boundary value problem completely specifled by choos- 
ing /i, q, fl, and ci. 

This solution clearly possesses quadrupolar symmetry, 
as deflned above. We have checked this solution numeri- 
cally for many values of the parameters. In addition, our 
numerical analysis did not flnd any other classes of solution. 
Thus, it appears that the analytic solution given above is the 
unique solution to our boundary value problem. We suspect 
that other classes of solutions, including singular solutions 
and solutions with dipolar symmetry, may exist if different 
boundary conditions are allowed. 

Not all of the solutions given by equations p9| to p7| 
correspond to protostellar collapse and outfiow solutions. 
(From equation ^ we note that outfiow solutions require 
Cl < and C2 > |ci|, so that material falls in near the 
equatorial plane and is ejected near the symmetry axis. With 
the help of equation ^ this restriction implies that ci < 
— 1. The opposite conditions provide us with solutions in 
which the fiow is reversed, with infall near the polar axis and 
expulsion near the equatorial plane. The maximum value 
of the self-similar dynamical pressure is cs + c^. Thus, we 
require cs -l-ce < 0, which keeps both the dynamical pressure 
and the out ward pressure gradient negative, as discussed in 
Section 



2.1 



There are, in fact, solutions with p > so that 
the pressure increases with radius. Such solutions might be 
physically interpreted as anisotropic collapse solutions that 
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are driven primarily by a sudden increase in the external 
pressure. This sort of collapse might be reasonable in regions 
where star formation is initiated by a strong external shock. 
However, we restrict ourselves to less exotic cases in this 
paper, where the pressure decreases radially outward. 

4.2 Characteristics of the Model and Exploration 
of the Parameter Space 

In this section, we discuss the general properties of the ana- 
lytical solution given by equations ^ to We supplement 
this discussion by a rather complete exploration of the pa- 
rameter space, in which we allow the parameters q, Q, and 
ci to vary within pre-defined ranges given by 



< M < 10 
-10^ < n < 10^ 



-10'' 
-10^ 



< g < 10'' 

< Ci < -1. 



(48) 



We sampled these ranges normally in each parameter pi 
when \pi\ < 1 and log- normally when \pi\ > 1. We find 
acceptable solutions right up to the edges of our parameter 
space given by equation Wq. In principle, we could expand 
the region of our parameter space search even further than 
the ranges given above. However, our choice to limit our pa- 
rameter space is consistent with the spirit of our self-similar 
analysis, where non-dimensional parameters are not usually 
expected to be too many orders of magnitude greater or less 
than unity. 

The Alfven singularity in our self-similar equations oc- 
curs when mA,e ~ 1, where 



rriAa 



(49) 



Note that the usual Alfvenic singularity is modified by 
our self-similarity ansatz so that only the components of 
the velocity and magnetic field are involved (eg. Tsinganos 
et al. 1996, LHF). This can be simplified with the help of 
equation pM 



2 



„2 • 



(50) 



The Alfvenic mach number has a constant value for each 
solution, which does not vary with either radius or angle. 
Both super-Alfvenic and sub-Alfvenic models are allowed, 
although there are no models where material passes through 
an Alfven point. Not, however, that the special solution 
where mA,e ~ 1 everywhere is allowed. Somewhat arbitrar- 
ily, we eliminate any solution for which mA,B < 0.1, which 
corresponds to magnetic energy densities in be greater than 
50 times the kinetic energy density in We- This provides 
us with a convenient way to exclude strongly sub-Alfvenic 
models with unrealistically high magnetic field strengths. 

There are several important conclusions that can be 
drawn directly from equations ^ t o ^From equation |39[ 
and with the help of equation^44|, the ratio of the outfiow 
velocity to the infall velocity is given by 

2 



WriO) 



Cl 



(51) 



Wri7T/2) 

The ratio WriO)/\Wr{ir /2)\ < 2 for aU models, since ci is 
negative. This upper limit on the net acceleration of the gas 
means that the axial outfiow velocity is at most twice the 
equatorial infiow velocity on a sphere of some given radius, 




Figure 1. We show an overlay of 6 typical solutions. Each model 
is colour coded for easy comparison between panels. 



at some instant of time. However, it is important to note 
that the velocity increases linearly with radius in our model, 
which is very different from the model presented in F Hl, 
where the velocity decreases with radius as r~^^^. Section 4.4 
shows that any given fluid particle escapes to radial infinity 
along the outflow axis in the finite time required for the 
central singularity to develop (at t = 0). Thus, large outfiow 
velocities are attained at some distance from the origin, near 
the time when the singularity develops. Realistically, the 
ejected material must encounter a shock with the external 
region at some finite radius, where the acceleration would 
presumably end. Nevertheless, the possibility of a strong jet- 
like outfiow remains, provided that the external medium is 
encountered at sufficient distance from the origin. 

The early stage of outfiow modelled in this paper is 
dominated by a quadrupolar circulation pattern with an 
outflow along the symmetry axis. It is very likely that the 
outflow would become even more vigorous at a later stage 
of evolution, in which energy is injected directly into the 
gas in the form of heat from the protostar. We have previ- 
ously included such heating in a steady state version of our 
quadrupolar flow model (FHl, LHF) and found that much 
higher velocities can indeed be achieved. Thus, the present 
model might represent the immediate predecessor in the evo- 
lutionary sequence leading up to our steady-state models. 

Figure |l] shows 6 of our solutions, chosen at random, 
overlaid on a single plot to illustrate the large dynamic range 
of self-similar velocities, magnetic field strengths, and pres- 
sures allowed. The curves are colour-coded for easy compar- 
ison between panels. The large allowed range of initial con- 
ditions and protostellar properties refiects the robustness of 
our model and of quadrupolar fiow in general. 

Now we examine the ratio of the radial magnetic field 
at the polar axis to the radial field at the equator. By com- 
bining equations ^ |39| and it is easy to show that 

Brio) 



Br{7T/2) 



= -2 



(52) 



for all solutions. This surprising result indicates that the fiow 
is structured such that the radial component of the magnetic 
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field is always compressed to precisely the same extent. It is 
interesting that we first noticed this result in outflow solu- 
tions obtained numerically using a shooting code. This was 
one of the key results that suggested the existence of a gen- 
eral analytic solution to our boundary value problem, and 
led us to the solution presented in equations to ^ 

A careful examination of Figure ^ reveals that most of 
our models have a low pressure region near the symmetry 
axis, but some demonstrate the opposite behaviour. This 
effect can be explained by consulting equations 
and Note that ce is the self-similar pressure difference 
between the midplane and the polar axis: 

C6 =p(7r/2)-p(0). (53) 

Note also that the term [aq + Q.) in equation ^ is the maxi- 
mum value of the self-similar toroidal field b^. Thus, a strong 
toroidal field tends to make cs negative, so that the pressure 
is highest at the polar axis. Clearly, this compression is due 
to the radial pinch of the toroidal field toward the polar 
axis (e.g. Fiege & Pudritz 2000a,b). On the other hand, the 
constant C4 in the second term is equal to the maximum 
rotational velocity, which acts in the opposite sense to push 
material away from the axis. This centrifugal barrier pro- 
duces a low pressure region near the axis when the toroidal 
field is not sufficiently strong to counter its effect. We refer 
to these models as "tornado" solutions, while we often refer 
to dominated models as magnetically pinched. A useful 
parameter to investigate these two types of behaviour is the 
ratio of the dynamical pressure at the the symmetry axis 
divided by the pressure at the midplane: 



P(0) 



C5 



C5 + C6 



(54) 



Since p < 0, tornado-type solutions correspond to 
p(0)/p(7r/2) > 1, while magnetically pinched solutions have 
the opposite behaviour, with p(0)/p(7r/2) < 1. In Figure ^ 
we plot this quantity against the magnetic pressure ratio 
brido)^ /bpoi{Go)^ , where fer(^o) and bpoi{6o) are respectively 
the self-similar toroidal field and poloidal field (defined by 



Opol 



y^br + bg) components, evaluated at the angle 



y/-ci/c2 



(55) 



which divides the equatorial infall zone from the polar out- 
flow zone (See equation Since ci < —1, equation ^ im- 
plies that Oo has a maximum possible value of cos^^ -^/l/S ~ 
54.7°. Note that the transition between tornado and mag- 
netically pinched solutions occurs when 6^,mai/^r,mai ~ 1, 
as one would expect. It is clear from the flgure that most of 
the parameter space is flUed with solutions of the tornado 
type. Both super and sub-Alfvenic solutions are represented 
in this figure, with no apparent discontinuity in behaviour. 
Note the sharp boundaries that are visible in the distribu- 
tion of tornado type solutions. These are due to two limits 
that can be shown analytically from ^ with the help of 
equations H and |4^: p(0)/p(7r/2) ^ 4 as ci —00, and 
p{0)/p{n/2) — > 1 as /i — > 00. 

It is interesting to comment on the origin of the toroidal 
field component in our model. The integral given in equation 
^ shows that the toroidal magnetic field can be expressed 
as a the sum of two terms: 



10' 



10" 



rIO - 



10 
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"Tornado" Solutions 



Pinched Solutions 



10" 



10" 



10"" 
b fevb ,(ej 

((1^ 0' pol^ 0' 



10= 
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Figure 2. Scatter plot of the ratio of the pressure at the axis 
over the pressure at the midplane, as a function of 6^ ^ax/^r,"iax 
for 5000 models. p(0)/p{n/2) > 1 corresponds to "tornado" type 
solutions, while p(0)/p(0)/p(7r/2) < 1 models are pinched by the 
toroidal magnetic field. (Recall that p < - see Fig. 1). 



Wr + 2/3 



Br + nsine 



(56) 



The first term is associated with toroidal field generated by 
VK0 twisting poloidal field lines, while the second term can 
be identified with toroidal fiux loops simply advected by the 
gas motions. 



4.3 Streamlines and Fieldlines 

Figure ^ is a split-frame figure showing the instantaneous 
stream (right hand side) and field (left hand side) lines in 
the poloidal plane, overlaid with vectors whose lengths are 
proportional to the magnitude of the poloidal velocity or 
magnetic field at each point. The stream and field lines are 
respectively given by solutions to the differential equations 

(57) 



r de 

dr 
VdB 



Wg 
Br 



(58) 



It is clear from the figure that the stream and field lines 
in the poloidal plane are not parallel. The stream and field 
lines are superimposed over a colour representation of the 
pressure on the left, and the toroidal field strength on the 
right. 

The path traced out by a fluid element in a time de- 
pendent model may not coincide with the instantaneous 
streamlines because the streamlines themselves evolve in 
time. However, an interesting property of our self-similar 
model is that a given fluid element followed in time actually 
obeys equation and hence moves only along the instanta- 
neous streamlines. This behaviour arises because time enters 
into each velocity component with the same power in equa- 
tion and the components of W do not explicitly depend 
on time. Thus the streamlines in the left panel of Figure 
^ can also be interpreted as the paths of individual fluid 
elements. 

Figure ^ shows a three-dimensional representation of a 



Protostellar Evolution during Time Dependent, Anisotropic Collapse 9 



N ol 




Figure 3. We show poloidal field fines (left) and streamlines 
(right) for a typical solution. These are overlaid on a colour plot 
of the toroidal field "pressure" b^/8iT on the left, and the gas 
pressure on the right. This particular model is a "tornado" type 
solution. The parameters of the model shown are as follows: /i = 
1.77, q = -6.11, n = -14.56, ci = -4.33. 




Figure 4. We show a three-dimensional view of several stream- 
linos (blue) and field lines (green). Note that there is considerable 
rotation in fiow, and that the streamline and field lines are not 
parallel, indicating the presence of a Poynting flux. The model 
shown is the same as the one in Figure ^. 



several stream and field lines equally spaced in angle and re- 
flected about the midplane. We have selected a model with 
a significant amount of rotation for this figure, whose pa- 
rameters given in the caption of figure ^ . The stream and 
field lines wraps around the origin as pressure gradients de- 
flect the material out along the symmetry axis. Note the 
similarity to the corresponding plot (Figure 5) in FHl. 



4.4 The Trajectories of Fluid Elements 

Some additional insight into the flow structure can be ob- 
tained by directly integrating the trajectories followed by 
fluid elements. It is trivial to integrate equation ^ with the 
help of equation ES, to show that 



tan 9 



-2C3 



(59) 



tan 9o 

where 9o is the poloidal angle of a fluid particle at some 
initial time to < 0. Combining equations ^ and ^ we find 
that 



2 + 3ci 
C3 = < 



-1/2 



(60) 



since ci < — 1 for all outflow solutions (Refer to the discus- 
sion following Thus, all fluid elements with 9a < 7r/2 
go to the symmetry axis 9 = when the central density be- 
comes singular at i = 0. Fluid elements on the asymptotic 
streamline that starts at 9o = n/2 go to the origin at t = 0. 
Note that an initially spherical shell of fluid is transformed 
by the flow into an thin "needle-like" distribution of material 
along the symmetry axis. 

It is also possible to obtain an exact expression for the 
radial part of the trajectory by directly integrating equation 
[s^ , using the forms given in equations |9| and We find 
that 



ro Vt 



tan 61 x^i/f^'^a) r 



/ sin 6* \ / cos 6*0 + 1 \ 
V sin 6*0 / V cos e + 1 y 



C2/{2C3) 



.(61) 



, tan 9o 

Therefore, r oc 6l''=i+'^2>'''^'=3) ^gg^j. ^^le time t = 0, when 
^ ^ 0. The discussion following equation ^ demonstrated 
that Ci + C2 > for all outflow solutions, and equation ^ 
indicates that C3 < 1. Thus, we find the intriguing result 
that material escapes to infinity for all outfiow solutions in 
the finite time required for the the central singularity to 
develop. This confirms that our model can indeed generate 
strong, jet-like outflows, provided that they can escape to 
sufficiently large distances before encountering a shock with 
the external medium. 



5 DISCUSSION 

In this section, we deal mainly with the question of where 
the preceding solution may be expected to apply. The reader 
may find it helpful to refer to Figure ^ throughout this dis- 
cussion, where we provide a sketch of the magnetic field 
configuration that we describe below. 

We proceed first with a qualitative discussion. For sim- 
plicity, we assume that the magnetic and angular momentum 
axes of a cold cloud core coincide. The collapse of a molecu- 
lar core may be due to the inside-out collapse of a SIS (Shu 
1977; GS), albeit one that is threaded by a magnetic field. 
Alternatively, it may be due to the loss of magnetic support 
due to ambipolar diffusion, or due to the damping of Alfven 
waves or MHD turbulence. The collapse might even result 
from the arrival of a sudden shock or torsional Alfven wave 
(Habe et al. 1991). Whatever the cause of the collapse, it 
is likely that a nearly radial magnetic field develops during 
collapse, as the field lines are dragged in with the collaps- 
ing gas. These radial field lines will be oppositely directed 
across the equatorial plane. It is in this plane (presumably 
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near the Keplerian point where the gas velocity may be rel- 
atively slow compared to the magnetic diffusion speed) that 
one can expect the central magnetic field to disconnect from 
the external field through the formation of an X type neu- 
tral point near the equator (e.g. Shu et al. 1994). The field 
lines should remain roughly radial outside of the equato- 
rial X points, especially near the midplane, as illustrated in 
Figure |5[ 

We must, however, also consider the pressure driven 
bounce of material in opposite directions above and be- 
low the equator along the axis. This rebounding material 
forms a quadrupolar velocity field (FHI; Tomisaka,1998; 
LHF), which distorts the radial field locally into rising axi- 
symmetric poloidal arches if the fiow is super- Alfvenic. As 
usual there will also be a toroidal component of the mag- 
netic field in these arches. These arches will also advect 
the toroidal component of the field with the rising gas. 
The key point is that these arches should develop first at 
small radii, where we expect the most vigorous bounce. Ul- 
timately, we expect this fiow pattern to lead to the forma- 
tion of quadrupolar velocity and magnetic fields consistent 
with the field structure that we have assumed in this series 
of papers. However, a full quantitative analysis remains to 
be done. The arches may themselves reconnect sporadically 
(through the formation of toroidal X lines), if the gas has 
finite resistivity so that ideal MHD does not strictly apply. 
This would have the effect of producing magnetic toroids 
looped around the axis, which would rise with the mean flow 
velocity. After each such disruption the persistent quadrupo- 
lar flow should reestablish the quadrupolar fleld. These mag- 
netic ejections may be associated with the production of 
energetic particles and hence "flaring" (e.g. Montmerle et 
al. 2000 for the detection of significant X-ray emission from 
class I protostars), and are probably of interest in their own 
right as features in the bipolar outfiow. 

The solution presented in this paper applies at smaller 
radii than the radius of the outer X point. As we have dis- 
cussed above, the bounce at this outer point is created by a 
mixture of rotational-magnetic support and radiative heat- 
ing (e.g. LHF), whereas the bounce at the inner X point 
is caused by the forming hydrostatic core. At this inner X 
point, we expect a transition to a dipolar magnetic field (al- 
though a quadrupolar circulation may still be present). Note 
that our equations may admit dipolar solutions, which could 
apply to the hydrostatic core at very small radii. However, 
we have not yet found any such solutions. Thus the settling 
solution that we have presented here is expected to lie be- 
tween inner and outer X points. Near the equator the field 
lines are actually wrapped around a magnetic O point. This 
equatorial XOX configuration may be a quite general mag- 
netic structure necessary to the process of star formation. It 
links the outer "Keplerian" disc to a more slowly growing 
inner hydrostatic core. The core probably grows to become 
the star at the expense of the 0-type envelope. 

One can also get some feeling for the applicability of 
our solution by looking at some numerical relation s. We 
can rewrite equation ( pi| ) to give t — {^/Gp) which 
shows that gives the time before the final singularity 
in units of the (approximate) free-fall time. Let us sup- 
pose that the region in which we are interested extends in- 
ward from about 1 au to the boundary of the hydrostatic 
core. If this region has a mass m in solar mass units then 




Figure 5. We sketch the possible regions where quadrupolar 
stream and field lines may arise during the formation of a star. 
The solution presented in this paper is concerned with the in- 
ner region while that presented in our earlier papers concerns the 
outer region. 



p « 2mx 10~'^(ji cm~^. Using m = O.IMq for a low mass star, 
we arrive &t t ~ ^JJi years. /,From equations ( |29[ ) and the 
solution (^) , we see that the inward equatorial velocity at 1 
au is about 4,ci/ y/p, km . Therefore if ci = —CyJJl where 
c is of order one but is sufflciently large to make p < (see 
equation then the bipolar velocity approaches 8c km 
and the equatorial velocity is comparable to the free fall ve- 
locity. 

It is also of interest to consider the magnitude of the 
r tiag netic field in this illustration. At 1 au the equation 
( pO|) plus the solution and our preceding assignments of the 
constants gives the equatorial radial field to be of order 
100 gauss. Of course this ignores the total magnetic fleld 
which is a bit more difficult to estimate because of the (j> 
component and its relation to the rotation of the system. 
Note that these numbers are provided for illustrative pur- 
poses only. Our models admits a large range of solutions 
that can describe a variety of protostellar objects. 



6 CONCLUSIONS AND SUMMARY 

We have developed a new model for the time dependent and 
anisotropic collapse that occurs within the inner regions of a 
star-forming molecular cloud core. By including the effects 
of self- gravity and MHD, our model provides a reasonably 
complete description of the dynamics on all scales between 
the inner hydrostatic core and an outer X point. Our previ- 
ous steady-state version of the model is expected to apply 
external to the collapsing region modelled here, and possibly 
at later times. Remarkably, the collapse model presented in 
this paper admits an exact and completely analytic solution. 
We note that there are few other analytic solutions of this 
complexity in all of MHD. We summarize our solution here, 
for the reader's convenience: 

Vr = — j(ci + C2 COS^ S) (62) 

ve = -^[C3sin(2e)] (63) 

T 

= --{c^smO) (64) 
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2 2 

1 

r 



Gt2 
r 



Gt^ 



q{ci + C2 cos^ ^ ^ 3^ 
qC3 sin(26') 



VGf- 
2 



■(gc4 + f2) sin 6* 



(65) 
(66) 
(67) 
(68) 
(69) 

(70) 



where the four free parameters of the model are fi, q, Q, and 
ci. The remaining constants in these expressions are given 
by equations ^ to Note that the time t starts large and 
negative, and the model evolves until a central singularity 
forms at t = 0. 

The main point of this work is to demonstrate that in- 
fall and outflow can coexist and arise naturally from our 
self-similar equations, with few additional assumptions. The 
outflow that arises is surprisingly vigorous, despite the lack 
of explicit internal heating in our model. We flnd that the 
outflow velocity increases linearly with radius, and that ma- 
terial in the outflow escapes to radial infinity in the finite 
time required for the central singularity to develop. More 
realistically, we expect the outflow to interact with the sur- 
rounding gas, outside of the self-similar region. This would 
limit the velocity of the outflow, but high velocities could 
still be obtained if the shock occurs much further out than 
the radius from which the material originated. 

Our model applies only at intermediate scales between 
two X type magnetic neutral points. We expect a transi- 
tion to a dipolar field internal to the inner X point, where 
the growing hydrostatic protostellar core resides. This re- 
gion could, in principle, be modelled using the same set of 
equations we have used in the protostellar collapse region 
modelled here, with boundary conditions consistent with a 
dipolar field. The region exterior to the outer X point is 
probably characterized by longer dynamical timescales than 
the protostellar collapse region. Once the overall quadrupo- 
lar fiow structure is established by the collapse, the outer 
region might be most appropriately modelled by the steady- 
state (or perhaps quasi-steady) version of this model, which 
we have previously discussed in depth (FHl, LHF). Future 
work on this class of models may proceed in two alternate 
directions. We may try to join together the three regions dis- 
cussed above, in piecewise fashion, to provide a more com- 
plete description of the simultaneous infall and outfiow pro- 
cesses. Alternatively, we might turn to the more ambitious 
problem of solving the most general sel f-similar partial dif- 
ferential equations given in Section 2.1 (equations ^ to [l8| ) 
over the entire range of scales. 
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